function [f,h,c] = WLD(image, T, M, S, ALPHA)
% WLD returns the local texture pattern of an image.

%% Determine the dimensions of the input image.
if size(image, 3) ~= 1
    image = rgb2gray(image);
end
d_image=double(image);
[ysize, xsize] = size(image);

%% Block size, each WLD code is computed within a block of size bsizey*bsizex
bsizey=3;
bsizex=3;

%% Minimum allowed size for the input image depends
% on the radius of the used LBP operator.
if(xsize < bsizex || ysize < bsizey)
  error('Too small input image. Should be at least (2*radius+1) x (2*radius+1)');
end

if nargin < 2
    T = 12;
end

if nargin < 3
    M = 4; %number of segment
end

if nargin < 4
    S = 4; %number of bin in each segment
end

if nargin < 5
    ALPHA = 6; % like a lens to magnify or shrink the difference between neighbors
end

%% filters
%    1  2  3  4  5  6  7  8  9
f00=[1, 1, 1; 1,-8, 1; 1, 1, 1];
f01=[0, 0, 0; 0, 1, 0; 0, 0, 0];
f10=[0,-1, 0; 0, 0, 0; 0, 1, 0];
f11=[0, 0, 0; 1, 0,-1; 0, 0, 0];

%% Compute the WLD code per pixel
%step 1 compute differential excitation
v00 = conv2(d_image,f00,'valid'); %convolve with f00
v01 = conv2(d_image,f01,'valid'); %convolve with f00
d_differential_excitation = atan(v00./v01); %perform the tangent scaling

%step 2 compute gradient orientation
v10 = conv2(d_image, f10, 'valid');
v11 = conv2(d_image, f11, 'valid');
c.c1 = nnz(v10 == 0);
c.c2 = nnz(v11 == 0);
c.c_and = nnz((v10 == 0) & (v11 == 0));
c.c_or = nnz((v10 == 0) | (v11 == 0));
c.C_AND = ((v10 == 0) & (v11 == 0));
d_gradient_orientation = atan2(v10, v11) + pi;
%image_GO range [0 .. 2*pi] -> [0 .. T-1]
image_GO=floor(d_gradient_orientation/(2*pi)*T);
image_GO(image_GO == T) = 0; 
h = hist(image_GO(:), T);

%% compute histogram
image_DE=wcodemat(d_differential_excitation, M*S, 'mat', 0);
f = hist3([image_DE(:),image_GO(:)], [M*S, T]) / numel(image_DE);
if nargout == 0
    hist3([image_DE(:),image_GO(:)], [M*S, T]);
    axis on
end
f = f(:)';
end %WLD_new
